/////////////////////////////////////////////////////////////////////
//
// \author Jeongnim Kim  jnkim@ncsa.uiuc.edu
//
// Measuring abstraction penalty using TinyVector and Vector based on PETE
// src/PETE      (copy of PETE library http://www.acl.lanl.gov)
// src/OhmmsPETE (Ohmms implementation of a fixed vector and an array)
//////////////////////////////////////////////////////////////////////
#include <iostream>
#include <vector>
#include <iomanip>

#if defined(ADD_)
#define addmult addmult_
#endif

#include "OhmmsPETE/TinyVector.h"
#include "OhmmsPETE/OhmmsVector.h"
#include "Utilities/RandomGenerator.h"
#include "Utilities/Timer.h"
#include "Utilities/OhmmsInfo.h"
using namespace APPNAMESPACE;
using namespace std;

int print_timing(const std::vector<double>& timingdata, int n, int niter, const std::string& what);

template<class C, class RNG>
class RandomVector{ };

template<class T, unsigned D, class RNG>
class 
RandomVector<TinyVector<T,D>, RNG> {

  public:
  typedef TinyVector<T,D> Return_t;
  RandomVector(RNG& gen): RG(gen){ }
  inline Return_t operator()() {
    Return_t res(0);
    for(int i=0; i<D; i++) res[i] = 2*RG()-1;
    return res;
  }
  RNG& RG;
};

// specialized for TinyVector<T,D>
template<class RNG>
struct RandomVector<TinyVector<double,3>, RNG> {
  typedef TinyVector<double,3> Return_t;
  RandomVector(RNG& gen): RG(gen){ }

  inline Return_t operator()() {
    return Return_t(2*RG()-1, 2*RG()-1, 2*RG()-1);
  }
  RNG& RG;
};

extern "C" { // fortran call
  void addmult(double* a,const double* b,const double* c, const double *d,
	       int* ntot);
}

inline void 
caddmult(double* restrict a, const double* restrict b, 
	 const double* restrict c, const double* restrict d, int ntot) 
{
  for(int i=0; i<ntot; i++) {
     a[i] = b[i] + 1.3*c[i] + 1.7*d[i];
  } 
}

int main(int argc, char **argv) {

  OhmmsInfo ohmms(argc,argv);
  std::string help("pete/ohmmspete");

  typedef TinyVector<double,QMC_DIM> Pos_t;
  typedef Vector<Pos_t> PosArray_t;
  typedef Vector<double> ValArray_t;

  int n=10;
  int i=0;
  int niter=5;
  while(i<argc) {
    std::string a(argv[i]);
    if(a == "-n"){
      n = atoi(argv[++i]);
    } else if(a == "-i") {
      niter = atoi(argv[++i]);
    } else if(a == "-c") {
      help = argv[++i];
    }
    i++;
  }

  if(argc<2) {
    cout << endl << "Using default size = " << n << endl;
    niter = 1000;
    cout << endl << "usage: -n size -i iteration -c keyword" << endl;
  }
  cout << "<!-- " << help << endl;
  cout << "Starting vector operations:" << n << " " << QMC_DIM <<"-dimensional elements" << endl;
  Random.init(0,1,-1); 
  RandomVector<Pos_t, RandomGenerator_t> rvgen(Random);

  double norm = 1.0/static_cast<double>(niter);
  Timer myclock;
  vector<double> timingdata(10);

  myclock.restart();
  PosArray_t A(n), B(n), C(n), D(n);
  timingdata[0] = myclock.elapsed();

  myclock.restart();
  A = 0.0;
  for(int i=0; i<n; i++) {
    B[i] = rvgen();
  }
  for(int i=0; i<n; i++) {
    C[i] = rvgen();
  }
  D = 1.1*B;

  timingdata[1] = myclock.elapsed();

  myclock.restart();
  for(int iter=0; iter<niter; iter++) {
    A = B + 0.13*C + 0.017*D;
  }
  timingdata[2] = myclock.elapsed()*norm;

  myclock.restart();
  for(int iter=0; iter<niter; iter++) {
    for(int i=0; i<n; i++)
      A[i] = B[i] + 0.13*C[i] + 0.017*D[i];
  }
  timingdata[3] = myclock.elapsed()*norm;
 
  if(n<=100) {
    for(i=0; i<n; i++)   {
      cout << A[i] << " " << B[i] << " " << C[i] <<  " " << D[i] << endl;
    }
  }

  ValArray_t a(QMC_DIM*n), b(QMC_DIM*n), c(QMC_DIM*n), d(QMC_DIM*n);
  int ii=0;
  for(int i=0; i<n; i++) {
    for(int j=0; j<QMC_DIM; j++)
      b[ii++] = B[i][j];
  }

  ii=0;
  for(int i=0; i<n; i++) {
    for(int j=0; j<QMC_DIM; j++)
      c[ii++] = C[i][j];
  }

  ii=0;
  for(int i=0; i<n; i++) {
    for(int j=0; j<QMC_DIM; j++)
      d[ii++] = D[i][j];
  }

  myclock.restart();
  for(int iter=0; iter<niter; iter++) {
    a = b + 0.13*c + 0.017*d;
  }
  timingdata[4] = myclock.elapsed()*norm;

  int ntot = QMC_DIM*n;
  myclock.restart();
  for(int iter=0; iter<niter; iter++) {
    for(int i=0; i<ntot; i++) a[i] = b[i] + 0.13*c[i] + 0.017*d[i];
  }
  timingdata[5] = myclock.elapsed()*norm;

  myclock.restart();
  for(int iter=0; iter<niter; iter++) {
    caddmult(&(a[0]), &(b[0]), &(c[0]), &(d[0]), ntot);
  }
  timingdata[6] = myclock.elapsed()*norm;

  myclock.restart();
  for(int iter=0; iter<niter; iter++) {
    addmult(&(a[0]), &(b[0]), &(c[0]), &(d[0]), &ntot);
  }
  timingdata[7] = myclock.elapsed()*norm;


  cout << "--> " << endl;
  print_timing(timingdata,n,niter,help);
}

/***************************************************************************
 * $RCSfile: vector.cpp,v $   $Author: jnkim $
 * $Revision: 1.5 $   $Date: 2005/12/11 15:28:52 $
 * $Id: vector.cpp,v 1.5 2005/12/11 15:28:52 jnkim Exp $
 ************************************************************************/
